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1.0  Project  Goals 


Our  research  investigated  the  feasibility  of  performing  image  enhancement,  based  upon 
diffusion  processing,  on  a  purely  quantum  or  hybrid  classical-quantum  computer.  In  the 
case  of  hybrid  computing  using  quantum  lattice  gas  algorithms,  the  approach  provides  a 
technological  road  map  for  reaching  quantum-only  computing.  A  hybrid  classical- 
quantum  processor  could  provide  an  interim  implementation  wherein  additional  qubits 
could  be  inserted  into  the  nodes  of  the  lattice  and  other  improvements  such  as  error 
correction  added.  In  this  way,  hybrid  classical-quantum  computing  could  serve  as  a 
testbed  for  a  given  hardware  approach. 

In  recent  years,  image  processing  based  upon  solving  partial  differential  equations 
(PDEs)  [perona,  lions,  sapiro,  osher]  has  become  more  prevalent  and  developed.  These 
techniques  include  digital  image  enhancement  and  multi-scale  image  representation  by 
way  of  nonlinear  PDEs.  We  are  motivated  to  examine  these  techniques  because  the 
fundamental  equation  of  non-relativistic  quantum  mechanics,  the  Schrodinger  equation, 
may  be  viewed  as  diffusion  (or  heat  conduction)  in  imaginary  time,  and  diffusion  in 
particular  provides  a  type  of  digital  image  filtering. 

Quantum  lattice  gas  algorithms  are  able  to  simulate  a  variety  of  PDEs  [yepez],  including 
the  Navier-Stokes,  Boltzmann,  and  Dirac  equations  [ibb].  In  our  research,  we 
investigated  which  diffusion  and  related  Schrodinger  equations  may  be  simulated  with  a 
quantum  lattice  gas  algorithm  (QLGA).  Additionally,  we  investigated  how  a  QLGA 
could  be  extended  to  perform  selective  image  smoothing  in  support  of  image 
enhancement. 


1 


2.0  Approach 

In  the  mathematical  sense,  the  partial  differential  equations  we  simulated  with  QLGAs 
are  evolution  equations  of  the  form 


ut  =  f(u,ux,uxx),  (2.1) 

where  the  subscripts  denote  partial  differentiation,  for  example,  ut  =  duldi.  Therefore,  we 
solved  equations  that  are  first  order  in  time  and  second  order  in  space.  The  PDE  is 
supplemented  by  an  initial  condition  u(x,t  =  0)  =  u(x,0)  and  boundary  conditions  at  the 
endpoints  of  the  x-interval  to  completely  specify  the  problem.  We  used  periodic 
boundary  conditions,  so  that  for  a  one-dimensional  problem  on  an  interval  of  length  L, 
u(0,t)  =  u(L,t).  For  image  processing  applications  in  particular,  the  dependent  variable  u 
corresponds  to  pixel  intensities.  For  heat  transfer  problems,  u  represents  the  temperature, 
while  for  fluid  flow  problems  u  corresponds  to  either  the  mass  density  or  flow  velocity. 

A  quantum  lattice  gas  algorithm  is  a  quantum  version  of  a  classical  lattice  gas,  which  in 
turn  is  an  extension  of  classical  cellular  automata  [yepez,  doolen].  In  place  of  the  binary 
lattice  variables  of  a  classical  lattice  gas,  the  quantum  version  has  a  local  Hilbert  space 
describing  the  quantum  bit  (qubits).  In  the  classical  case,  in  order  to  recover  the  proper 
macroscopic  dynamics  (and  thermodynamics);  it  is  important  to  ensure  that  the 
microscopic  dynamics  preserves  conservation  laws.  Similarly,  in  the  quantum  case,  the 
number  densities  of  qubits  must  be  preserved  so  that  the  interaction  operator,  called  the 
collision  operator,  must  be  unitary. 

The  operation  of  a  type-II  quantum  processor  includes  the  sequential  repetition  of  three 
main  steps  [berman,  yepez,  vahala,  love].  First,  initialization  creates  the  quantum- 
mechanical  initial  state  that  corresponds  to  the  initial  probability  distribution  for  a  partial 
differential  equation  to  be  solved.  Secondly,  a  unitary  transformation  (collision  operator) 
is  applied  in  parallel  to  all  the  local  Hilbert  spaces  in  the  lattice.  Fastly,  in  the 
measurement  step,  the  quantum  states  of  all  the  nodes  are  read  out.  These  results  are 
used  to  reinitialize  the  quantum  processor  in  the  state  which  corresponds  to  the  new 
probability  distribution. 

After  the  initialization  step,  a  quantum  lattice  gas  algorithm  performs  iterations  of  a 
collision  operator  and  a  streaming  operator.  The  latter  operator  shifts  the  state  of  a  qubit 
from  a  given  lattice  site  to  its  nearest  neighbors  in  the  lattice. 

It  turns  out  that  both  of  the  streaming  and  collision  operations  in  a  QFGA  may  be 
combined  into  a  single  discretized  lattice-Boltzmann  equation  (FBE) 

fa(x+slea,t+s2x)  -  fa(x,t)  =  Qa(x,t),  (2.2) 

where  Qa(x,t)  =  v|/((x,t)|UlnaU  -  na|\|/(x,t))  is  the  collision  term.  Here  ea  denotes  the  lattice 
direction  for  streaming  of  qubit  a,  1  is  the  lattice  spacing,  t  the  time  increment,  and  na  the 


2 


multi-qubit  number  operator.  The  small  parameter  s  is  the  ratio  of  1  to  L,  the  size  of  the 
ID  system,  and  corresponds  physically  to  the  Knudsen  number  of  fluid  dynamics.  The  fa 
equation  above  is  an  analog  of  the  discretized  classical  Boltzmann  equation  for  a  particle 
distribution  function  f(x,v,t). 

The  collision  term  vanishes  when  U|\|/(x,t)>  =  |\|/(x,t)).  This  gives  the  equilibrium 
distribution  faeq. 

The  collision  operator  U  is  further  constrained  by  imposing  the  analog  of  mass  and 
momentum  conservation.  These  conditions  also  constrain  the  collision  terms  Qa. 

For  the  mass  density  p  =  Za=iB  mafa(x,t)  to  be  constant,  we  see  by  multiplying  the  LBE  by 
ma  and  summing  over  the  B  qubits  that 

La=!B  maQa(x,t)  =  0.  (2.3) 

Similarly,  multiplying  the  LBE  by  cmaeai,  where  c  is  a  constant  speed,  and  summing 
gives 


cZa=iB  maeaiQa(x,t)  =  0.  (2.4) 

This  represents  momentum  conservation  in  the  ith  direction. 

Parallel  to  the  mass  density  p  =  Za=iB  mafa  is  the  momentum  density 

pv;  =  cZa=iB  maeaifa(x,t).  (2.5) 

In  enforcing  mass  and  momentum  conservation  we  are  effectively  taking  moments  of  the 
lattice  Boltzmann  equation  to  derive  macroscopic  quantities  for  the  lattice  gas  system.  In 
quantum  mechanical  (QM)  terms,  the  conserved  quantities  have  corresponding  operators 
that  commute  with  the  collision  operator  U.  For  instance,  the  mass  operator  is  given  by 
Qo  =  Za=iB  mana.  Then  Za=iB  ma  Qa(x,t)  =  \j/((x,t)|UtQ0U  -  Q0|\|/(x,t)>  =  0,  as  previously 
described  for  the  equilibrium  condition. 

On  the  theoretical  side  of  our  efforts,  we  made  advances  in  both  deriving  finite  difference 
approximations  for  n-dimensional  linear  diffusion  and  in  the  so-called  Chapman-Enskog 
analysis  of  simpler  QLGAs.  For  multidimensional  anisotropic  linear  diffusion  we  now 
have  a  QLGA  based  upon  two  qubits  per  node  with  theory  and  simulation  verifying  the 
form  of  the  diffusion  coefficients. 

We  have  made  progress  in  developing  the  Chapman-Enskog  multi-scale  expansion  for 
ID  QLGAs  with  two  qubits  per  node.  This  type  of  expansion  requires  a  detailed  study  of 
the  associated  equilibrium  state  of  the  system,  requires  the  enforcement  of  mass  and 
momentum  conservation,  and  the  pseudo  inversion  of  a  Jacobian  matrix.  The  Chapman- 
Enskog  expansion  is  carried  out  in  the  small  parameter  s  that  corresponds  to  the  classical 
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Knudsen  number.  For  our  QLGAs,  this  corresponds  to  the  ratio  of  the  lattice  spacing  1  to 
the  number  of  nodes. 

The  Jacobian  matrix  is  a  result  of  linearizing  the  right-side  of  the  lattice  Boltzmann 
equation  with  the  collision  operator.  It  is  in  general  a  rank  deficient  matrix,  reflecting  the 
imposition  of  conservation  rules. 

The  Chapman-Enskog  analysis  tends  to  be  rather  involved,  even  for  simpler  ID 
algorithms.  It  additionally  includes  algebraic  inversions  for  the  equilibrium  occupation 
probabilities.  A  later  section  of  this  report  details  the  Chapman-Enskog  analysis  for  a 
specific  ID  algorithm  with  two  qubits  per  node  and  corresponding  to  unit  mass  particles, 
ma  =  1,  for  a  =  1  and  2. 

Under  the  next  heading  of  this  report,  we  present  details  first  for  simulation  of  the  linear 
diffusion  equation  and  then  for  Burgers’  nonlinear  equation.  Further  details  may  be 
found  in  an  accompanying  manuscript  that  has  been  submitted  for  external  publication. 
The  discussion  for  the  linear  diffusion  equation  includes  the  instance  of  an  in-depth 
Chapman-Enskog  analysis.  We  then  proceed  to  a  description  of  QLGA  simulation  for  the 
linear  and  nonlinear  Schrodinger  equations. 
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3.0  Details  of  Technical  Accomplishments 
3.1  Linear  Diffusion 


Diffusion  processing  has  proved  very  useful  for  practical  image  enhancement,  wherein 
the  visual  quality  of  an  image  is  improved.  We  have  investigated  methods  of  carrying  out 
such  processing  in  a  combined  classical-quantum  computing  environment.  While 
quantum  Fourier  and  some  quantum  wavelet  transforms  are  known  [nielsen],  their  direct 
application  to  image  or  signal  processing  is  unclear.  Motivated  by  real-world 
applications,  we  offer  an  alternative  based  upon  the  simulation  of  PDEs. 

Key  to  the  efficient  simulation  of  PDEs  on  a  hybrid  (or  Type  II)  processor  is  the 
execution  of  a  QLGA.  In  a  type-II  architecture,  nodes  of  phase-coherent,  entangled 
qubits  are  linked  by  classical  communication  to  nearest  neighbors  in  a  discrete  lattice 
[yepez,  yepezOl].  Generally  this  type-II  approach  has  the  advantage  of  requiring  both 
less  spatial  and  temporal  entanglement  than  for  the  usual  (type-I)  quantum  computing 
methods.  Indeed,  if  the  coherence  time  of  the  qubits  should  be  sufficiently  long,  quantum 
error  correction  would  not  be  required.  However,  such  error  correction  could  be 
implemented  within  the  nodes  if  needed,  or  as  a  testbed  of  an  interim  quantum  computing 
(QC)  technology. 

A  QLGA  includes  the  sequential  repetition  of  three  main  steps  [yepez,  yepezOl,  berman, 
vahala,  yepez06].  First,  initialization  creates  the  quantum-mechanical  initial  state  that 
corresponds  to  the  initial  probability  distribution  for  a  partial  differential  equation  to  be 
solved.  Secondly,  a  unitary  transformation  is  applied  in  parallel  to  all  the  local  Hilbert 
spaces  in  the  lattice.  Lastly,  in  the  measurement  step,  the  quantum  states  of  all  the  nodes 
are  read  out.  These  results  are  used  to  reinitialize  the  QC  in  the  state  which  corresponds 
to  the  new  probability  distribution. 

The  type-II  approach  offers  a  means  to  effectively  solve  complex  gas  and  fluid  dynamics 
problems  [yepez,  vahala].  While  QLGAs  have  been  shown  to  solve  a  number  of  PDEs, 
implementation  of  these  algorithms  has  been  mainly  done  for  one  spatial  dimension  only. 
Very  recently  multi-dimensional  simulation  for  the  Schrodinger  equation  has  been 
performed  [sakai]. 

Classical  lattice  Boltzmann  simulations  have  proved  highly  valuable  for  fluid 
phenomenon,  including  the  presence  of  complex  geometry.  Quantum  lattice  gas 
algorithms  may  be  written  in  an  analogous  discretized  lattice  Boltzmann  equation  form  as 
mentioned  in  Section  2.0. 

There  exist  quantum  lattice  gas  algorithms  that  do  not  correspond  to  the  average  over 
some  underlying  lattice-gas  model  in  the  Boltzmann  approximation  and  their  complexity 
analysis  remains  an  open  problem  [love].  The  algorithms  we  consider  do  have  a  classical 
counterpart  and  these  provide  an  efficient  implementation  on  parallel  architectures. 
Therefore  our  work  highly  suggests  the  usefulness  of  both  classical  and  hybrid  classical- 
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quantum  lattice  gases  for  image  enhancement.  Our  methods  are  well  suited  for 
implementation  on  parallel  computers,  and  beyond  this,  to  computational  grids  (e.g., 
[harting]). 

The  appeal  of  our  methods  for  classical  computers  is  especially  important  as  the 
computational  advantages  of  purely  quantum  processors  continue  to  be  elusive.  This  is 
illustrated  by  the  extremely  recent  and  significant  discovery  that  the  quantum  Fourier 
transfonn  may  be  efficiently  classically  simulated  [dorit] .  With  this  finding,  there 
appears  to  be  new  opportunities  for  novel  algorithms  for  classical  digital  and/or  digital- 
analog  computers.  In  at  least  some  sense  it  appears  that  the  complexity  characterization 
for  the  quantum  Fourier  transform  had  been  otherwise  anticipated  [alicki] . 

In  a  separate  preprint  [coffeycolburn],  we  present  simulation  examples  of  multi¬ 
dimensional  linear  diffusion  on  a  type-II  quantum  computer.  For  the  sake  of  definiteness, 
we  mainly  concentrate  on  two-dimensional  (2D)  diffusion  which  is  mapped  onto  a 
rectangular  lattice  of  nodes  in  a  type-II  QC.  However,  our  algorithm  carries  over  to 
higher  dimensions  and  we  briefly  describe  this.  We  may  stress  that  diffusion  processes 
are  common  throughout  engineering  and  scientific  fields.  Therefore,  diffusion  modeling 
is  valuable  amongst  a  very  large  number  of  disciplines,  extending  beyond  signal  and 
image  processing. 

Besides  the  sections  of  the  paper  [coffeycolburn]  dealing  with  2-,  3-  and  n-dimensional 
linear  diffusion,  we  present  an  important  innovation  for  image  applications.  Namely,  we 
demonstrate  constrained  linear  diffusion  that  provides  a  method  for  non-uniform  image 
smoothing.  This  extended  algorithm  fits  naturally  into  our  type-II  computing  approach, 
combining  a  constraint  condition  on  pixel  values  differences  with  the  measurement  step. 

We  have  implemented  a  variety  of  QLGAs  for  linear  diffusion  in  1-4  spatial  dimensions 
in  both  Matlab  and  Mathematica.  Of  these,  we  have  elsewhere  [coffeycolburn]  reported 
in  detail  the  Mathematica  implementations  for  2D  and  3D  simulation.  The  constrained 
2D  diffusion  algorithm  has  been  implemented  in  Mathematica. 

For  linear  diffusion  the  local  collision  operator  U  can  be  based  upon  the 
V SWAP  gate  [yepez,  bennan,  yepezOl], 


10  0  0 

o  ^(1-0  | (i  +  0  0 

o  |(i  +  o  | (i-0  o 

0  0  0  1 


(3.1) 


In  this  equation  there  is  a  U(2)  subblock  entangling  states  |01)  and  |10). 
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A  preliminary  implementation  of  a  type-II  QC  has  been  made  in  a  liquid-state  nuclear 
magnetic  resonance  (NMR)  system  [pravia],  These  experiments  used  a  chloroform 
(  CHCI3)  solution,  with  hydrogen  and  a  particular  carbon  nuclei  serving  as  two  qubits 
per  node.  Sixteen  nodes  of  a  one-dimensional  lattice  were  created  by  the  magnetic  field 
of  a  gradient  coil.  Specific  radio  frequency  (rf)  pulses  were  applied  for  the  unitary 
operations  on  qubits.  It  was  possible  to  carry  out  a  dozen  time  steps  of  the  algorithm  for 
the  linear  ID  diffusion  equation.  Improved  system  control  should  pennit  the  execution  of 
the  algorithm  for  longer  times  with  more  fidelity.  The  NMR  ensemble  computing 
approach  provides  a  fascinating  combination  of  quantum  computing  with  classical 
molecular  computing  and  it  has  been  estimated  that  implementations  with  30  or  more 
addressable  qubits  would  exceed  the  computational  capacity  of  any  classical  digital 
computer  (e.g.,  [love]). 

Many  further  avenues  for  physical  implementations  of  type-II  QCs  should  exist.  These 
include  other  spin-based  systems,  coupled  ion  or  neutral  atom  traps,  and  flux-,  charge-,  or 
phase-based  superconducting  qubit  systems  [herns].  An  attraction  of  the  type-II  approach 
is  the  possibility  of  realization  significantly  prior  to  large-scale  type-I  QCs. 


3.1.1  Constrained  Linear  Diffusion  Processing 

To  be  of  most  use  to  image  processing  tasks,  we  would  like  to  have  the  ability  to  perform 
selective  image  smoothing  and  linear  diffusion  alone  will  not  provide  this.  In  this 
subsection  we  describe  and  demonstrate  a  technique  that  gives  non-uniform  image 
smoothing  within  our  type-II  computing  approach.  We  combine  a  linear  diffusion  step 
with  a  constraint  condition  step.  Since  measurement  of  qubit  states  is  part  of  a  quantum 
lattice  gas  algorithm,  the  incorporation  of  the  constraint  step  fits  well  as  part  of  the 
overall  algorithm.  An  example  of  the  use  of  constrained  linear  diffusion  for  color  image 
dequantization  with  classical  computing  is  given  in  [keysers]. 

Suppose  that  our  original  image  is  Io.  We  now  perform  at  iteration  j,  j  =  1,2,. . .  on  image 

If 

(i)  2D  linear  diffusion,  giving  the  image  Ij',  followed  by 

(ii)  the  pixelwise  constraint 

Ij+i(x,y)  =  max[I0(x,y)  -  8,  min(I0(x,y)  +  8,  Ij'(x,y))].  (3.2) 

Here  8  is  a  processing  parameter.  It  can  be  chosen  according  to  the  distribution  of  pixel 
values  or  other  properties  of  an  image.  For  instance,  8  can  be  taken  as  proportional  to  the 
quantization  level  for  pixel  intensities.  For  gray  scale  images  with  pixel  values  in  the 
range  0  to  2  -1  we  could  choose  8cc  2'  . 

Note  that  the  measurement  and  constraint  step  in  the  above  algorithm  effectively  gives  us 
a  form  of  nonlinear  operation  within  the  processing.  This  is  an  advantage  of  the  type-II 
computing  approach.  With  a  pure  quantum  processor,  this  would  not  be  possible.  With 
the  constraint  in  place,  image  regions  with  similar  pixel  intensities  are  diffused  while 
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regions  with  large  differences  such  as  edges  remain  hardly  touched.  This  is  highly 
attractive  for  image  processing,  as  the  edge  information  is  preferentially  preserved. 

Trial  experience  with  this  constrained  linear  diffusion  algorithm  indicates  that  10-100 
iterations  are  effective.  This  fairly  modest  number  of  iterations  means  that  practical 
image  sizes  (of  a  few  hundred  by  a  few  hundred  pixels)  could  be  routinely  processed. 

Our  implementation  of  constrained  2D  diffusion  has  been  done  in  Mathematica  for 
images  with  8-bit  intensity  values.  Each  of  the  occupancy  probabilities  fj(x,y,0)  for  the 
two-qubit  algorithm  is  initialized  as  one -half  of  the  pixel  intensity  divided  by  255.  When 
displaying  intermediate  or  final  image  processing  results,  the  image  values  are  multiplied 
by  255  in  order  to  use  the  Mathematica  graphics  functions. 

In  order  to  have  a  benchmark  for  our  constrained  diffusion  algorithm,  we  have  taken  a 
bitmap  input  image  of  the  histogram-equalized  image  of  Figure  2b  of  [keysers] — see 
Figure  1.  This  input  image  contains  a  rectangular  region  of  smooth  grayscale  variation, 
rectangular  striped  regions,  and  a  sample  alphabetic  text  region.  Performing  linear 
diffusion  only,  the  image  is  quickly  uniformly  blurred.  This  is  demonstrated  in  Figure  2 
after  10  and  20  linear  diffusion  steps.  Shown  on  the  right  side  of  Figure  2  is  the 
difference  of  the  processed  image  and  the  original  image. 

Upon  the  test  image  of  Figure  1  we  have  additionally  run  50  steps  of  the  above 
constrained  algorithm  with  8  =  0.025.  The  left  column  of  output  images  at  successive 
iterations  in  Figure  3  shows  that  the  features  of  the  input  are  preserved  over  time.  In 
order  to  bring  out  the  pixel  changes  during  the  constrained  diffusion  processing,  in  the 
right  column  of  Figure  3  we  show  the  difference  between  the  original  image  and  the 
result  of  the  constrained  processing  at  the  corresponding  iteration.  According  to  the 
sequence  of  images  on  the  left  side,  we  see  that  the  desirable  property  of  intra-region 
smoothing  has  occurred.  For  instance,  in  the  largest  single  rectangular  region  in  the  top 
of  the  image,  the  original  grayscale  stepping  has  been  smoothed.  This  has  occurred  while 
the  striping  and  textual  infonnation  in  the  rest  of  the  image  has  been  maintained. 
Therefore,  this  algorithm  should  be  a  viable  candidate  for  image  pre-  or  post-processing 
or  image  enhancement  on  a  type-II  quantum  computer.  We  have  tried  the  constrained 
diffusion  algorithm  on  other  input  images  with  similar  evidence  of  intra-region 
smoothing. 


a  bcdefg  h  ij  klm  nopq  rstu  v  wxyz 


Figure  1.  Test  image  for  constrained  linear  diffusion  processing. 
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Figure  3.  Constrained  linear  diffusion  processing  example  after  1,  25  and  50  steps. 
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Figure  3  is  an  example  of  constrained  linear  diffusion  processing.  The  images  in  Figure  3 
show  the  results  after  1,  25,  and  50  steps.  With  a  constraint  step  (on  pixel  value 
differences)  included  in  the  overall  algorithm,  the  edge  infonnation  within  the  image  is 
maintained. 

By  way  of  summary  and  emphasis,  we  may  state  that  within  the  constrained  diffusion 
processing  approach,  we  effectively  achieve  nonlinear  filtering  of  a  signal  or  image, 
without  the  much  greater  complexity  of  solving  a  nonlinear  diffusion  equation.  The 
constrained  linear  diffusion  algorithm  contains  a  parameter  8  for  pixel  wise  comparison 
of  the  intensities  of  the  current  with  the  previous  and  original  images.  We  also 
investigated  other  variations  on  constrained  linear  diffusion,  including  the  choice  and 
possible  adaptivity  of  the  constraint  parameter  8. 


3.2  Chapman-Enskog  Analysis  for  Two-Qubit  Algorithm 

In  this  section  we  describe  the  Chapman-Enskog  analysis  for  a  QLGA  for  linear  diffusion 
with  two  qubits  per  node.  Suppose  that  the  lattice  Boltzmann  equation  takes  the  form 

fa(x+slea,t+s2T)  -  fa(x,t)  =  Qa(x,t),  (3.3) 

where  1  is  the  lattice  spacing,  ea  the  momentum  directions  for  the  particles  (ath  qubit),  t 
the  time  step,  and  Qa  the  collision  tenn.  For  this  relatively  simple  example,  a  =  1,  2  only 
and  ei  =  1,  e2  =  -1.  The  occupancy  probabilities  fa  give  the  probability  amplitudes  for 
each  of  the  two  qubits  as 


|qa(x,t)>  =  Vfa(x,t)}|  1)  +  V(l-fa(x,t))|0).  (3.4) 

The  (small)  expansion  parameter  s  ~  1/L,  where  L  gives  the  size  of  the  ID  system, 
corresponds  to  the  classical  Knudsen  number.  (1  serves  as  the  mean  free  path  between 
collisions  here.) 

In  this  case  the  collision  terms  may  be  written  as  Qi  =  -Q2  =  Ef  where 


Q  =  1/2[f2(l-f1)-f1(l-f2)]. 


(3.5) 


They  satisfy  Qi  +  F22  =  0,  corresponding  to  mass  density  conservation  for  unit  mass 
particles.  We  may  note  that  the  Jacobian  matrix  Jab  =  cQa  /of,, 


J  = 


jY-1 

2  V  1 


n 


(3.6) 


is  independent  of  the  occupational  probabilities  fa.  Equilibrium  occurs  for  Qa  =  0  or  fi  = 
f2=  d/2.  The  Jacobian  is  a  special  case  of  a  Hermitian  matrix.  Its  two  Eigen  values  must 
be  real  and  their  sum  ki  +  A,2  =  -1  =  Tr  J,  where  Tr  denotes  the  trace  operation. 
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Taylor  series  expanding  equation  (3.3)  gives 


2j2 


sleA+eAA+A 

a  dx  dt  2 


dx 2 


+  o(si)  =  na(L)  +  zl 


5A 


^+- 


(3.7) 


2 

where  ea  =  1  and  ...  represent  higher  order  terms.  For  the  Chapman-Enskog  expansion 
we  take 


fi=fi(0)+sfi(1)+s2fi(2)+  ...,  (3.8a) 

and 

f2=  f2(0)+  s  f2(1)+  s2  f2(2)+  ...,  (3.8b) 

where  f/0’  and  f2<0)  are  the  equilibrium  values  of  the  occupancy  probabilities.  Here  they 
are  particularly  simple  and  each  given  by  d/2.  The  first  term  of  the  right  side  of  equation 
(4)  vanishes  when  evaluated  at  equilibrium. 

In  this  two-qubit  example,  much  of  the  full  Chapman-Enskog  analysis  may  be  avoided. 
We  first  sum  equation  (4)  over  a  =1,  2  and  use  the  mass  density  p  =  fi+f2,  giving 


,  d  2  dp  s2/2  d2 p 

si  —  (f  -  f0)  +  szT  —  + - v-  =  0. 

dx  1  2  dt  2  dx2 


(3.9) 


Here  we  interchanged  the  double  sum  on  the  right  side  and  applied  conservation  of  mass. 

It  appears  that  the  first  term  on  the  left  side  of  equation  (3.9)  also  contributes  at  0(s“)  and 
we  proceed  to  work  out  the  details  for  s(fi(1)-  f2(1)).  Consider  the  O(s)  equations  from 
equation  (4).  Explicitly  these  read 


fe  A 1 _z 

‘  a*  A  dfh  b 


By  making  use  of  the  Jacobian  matrix  of  equation  (3.6)  these  equations  are: 

'(0)  j 


p\ 


dx  2 


and 


-/ 


dfi 


(0)  | 


.  (1)-/2(1)]- 
dx  2 


Summing  equations  (3.11a)  and  (3.1  lb)  gives: 

1  (d/dx)  (fi(O)-f2(0))  =  0. 


(3.10) 


(3.11a) 

(3.11b) 


(3.12) 
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This  is  expected  since  fi(0)  =  f2<0).  This  consistency  check  is  a  special  case  of  a  much 
more  general  solvability  condition  corresponding  to  general  equations  of  the  form  of 
equation  (3.10)  and  is  further  discussed  in  the  Appendix. 

The  solvability  condition  arises  due  to  the  fact  that  the  Jacobian  matrix  J  appearing  in 
equation  (3.10)  is  not  invertible.  Equation  (3.10)  may  also  be  written  in  matrix-vector 
form 


Qf (0) 

l  \e  ~ —  >=  J  I  /(1) 

a  OX 


>, 


(3.13) 


where:  If*1*)  = 


a) 


/i 

/2(1> 


.  The  Jacobian  has  two  Eigen  values  X\  =  0  and  k2  =  -1  with  (right) 


T  T  T  2 

Eigen  vectors  (1,1)  and  (-1,1)  respectively,  where  denotes  transposition.  Note  that  J" 
=  -  J  =  A.2  J.  So  multiplying  equation  (3.13)  on  the  left  by  J  gives 


df(0) 

U  I  ea  A  >=  hJ  I  f'”  >  . 

OX 


(3.14) 


We  take  this  equation  to  imply 


1  |ea  (3fa(0)/5x)>  =  A.2 1^1}), 


(3.15) 


omitting  discussion  of  the  pseudo  inverse  of  a  matrix.  Component  wise,  equation  (3.15) 
states  that 


1  (3fi(0)/3x)  =  -fi(1),  (3.16a) 

-1  (5f2(0)/Ox)  =  -f2(1),  (3.16b) 

i.e.,  the  first  order  corrections  fi(1)  and  f2(1)  are  directly  proportional  to  the  derivative  (the 
gradient  in  higher  dimensions)  of  the  corresponding  equilibrium  distribution. 

Now  consider  the  (first  nonzero)  contribution  at  0(s  )  from  the  first  term  on  the  left  side 
of  equation  (3.9)  that  we  may  identify  from: 


df 

sl'Lea  — -  =  sTLe 

Cl  Cl  ^  Cl  c 

ox 


'(2) 


K\SWL+S>£L+. 


dx 


dx 


dx 


(3.17a) 


e^aea 


dL 


(i) 


dx 


■  +  £ 


(2) 


dx 


■  +  . 


(3.17b) 
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d/r 

dx 


sir 

dx 


+  <sX 


5/a<2) 

dx 


+  ... 


,d2f! 

dx2 


,d2tt 

dx2 


-  +  0(e) 


+  0(e) 


(3.17c) 

(3.17d) 

(3.17e) 


Here  equation  (3.17b)  follows  from  equation  (3.17a)  by  momentum  conservation  (or, 
equivalently,  by  equation  (3.12)).  Equation  (3.17d)  follows  from  equation  (3.17c)  by 
applying  the  relations  (3.16).  We  conclude  that 


az.e^=d^-(fl-f,)  =  -sV 

dx  dx 


d2P 

dx2 


+  0(e) 


Then  equation  (3.9)  becomes  (at  0(s“)) 


r^+;2 


r  i  A 


dt 


--1 

v2  j 


d2p 


or 


dp 

~dt 


v2ry 


dt2 

d2p 


=  0, 


dt 


2  ? 


(3.18) 


(3.19a) 

(3.19b) 
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the  linear  diffusion  equation  with  diffusion  constant  D  =  1  /2t. 


3.2.1  Additional  Chapman-Enskog  Solvability  Discussion 

We  further  discuss  the  solvability  condition  for  the  first  order  equations  (3.10)  and  the 
null  Eigen  vectors  of  the  Jacobian  matrix.  We  assume  B  qubits. 

The  mass  density  conservation  condition  for  the  collision  terms  is  Za=iB  maOa(  { f, } )  =  0. 
Hence,  expanding  about  equilibrium,  fa  =  fa(0)+  s  fa(1),  we  have  the  condition  Za.b  ma  Jab 
fb<1)  =  0  where  fb(1)  is  arbitrary.  Thus  the  vector  (M|  =  (mi,m2,...,mB)  is  a  left  null  Eigen 
vector  of  J,  (M|J  =  0. 

Multiplying  equation  (3.10)  on  the  left  by  (M|  we  obtain  the  solvability  condition  lZa  ea 
ma  3fa(0V3x  =  0.  For  the  equal-mass,  two-qubit  example  above,  this  equation  is  trivially 
satisfied. 
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The  momentum  density  conservation  condition  for  the  collision  tenns  is  Za=iB  ma 
eaQa({fj})  =  0.  Expanding  about  equilibrium,  we  obtain  the  condition  Xa,b  ma  ea  Jab  fb<1>  = 
0  where  again  fb(1)  is  arbitrary.  This  equation  implies  (u|  =  (miei,...,mBeB)  is  also  a  left 
null  Eigen  vector  of  J.  This  condition  does  not  hold  for  the  one  dimensional  two-qubit 
example  above  where  the  average  momentum  is  zero. 

The  Jacobian  matrix  of  equation  (3.6)  is  special  in  another  sense.  It  is  an  example  of  a 
symmetric  circulant  matrix,  for  which  Jab  —  j|a-b|,  where  jc  is  a  vector  of  coefficients.  Since 
this  operator  is  real  and  symmetric,  it  has  only  real  Eigen  values.  For  the  QLGAs,  J  is 
generally  rank  deficient  and  so  has  at  least  one  zero  Eigen  value. 

Since  the  entries  of  Jab  depend  only  upon  the  difference  a  -  b,  J  is  also  a  Toeplitz  matrix. 
Very  many  results  are  known  on  the  Eigen  structure  of  Toeplitz  and  circulant  matrices. 
Typically  for  QLGAs,  the  coefficients  of  jc  are  highly  related  and  multiple  zero 
digenvalues  result. 

A  reference  for  the  Chapman-Enskog  expansion  for  classical  computational  fluid 
dynamics  is  the  book  Lattice-gas  Cellular  Automata  by  D.  H.  Rothman  and  S.  Zaleski, 
Cambridge  University  Press  (1997).  (See  especially  Ch.  15.) 


3.3  Burgers’  Equation 


In  this  section  subscripts  once  again  denote  partial  differentiation.  Burgers’  equation 


ut  +  uux  =  vuxx  (3.20) 

has  been  used  in  simplified  models  of  fluid  turbulence.  Below  we  recall  two  of  its  basic 
properties,  the  Hopf-Cole  transfonnation  to  the  linear  heat  equation,  and  the  energy 
integral. 

In  contrast  to  the  complex-valued  collision  operator  based  upon  the  square-root-of- 
SWAP  gate  used  for  the  linear  diffusion  equation,  for  Burgers’  equation  the  appropriate 
collision  operator  at  a  lattice  site  is  proportional  to  the  square-root-of-NOT  gate, 


10  0  0 
o  1/V2  1/V2  0 

0  -1/V2  1/V2  0 

0  0  0  1 


(3.21) 


This  matrix,  that  is  overall  unitary,  has  a  unitary  2x2  sub-block  mixing  only  the  states 
|01)  and  1 10)  and  thereby  preserves  particle  number  (probability).  This  matrix  can  be 
thought  of  as  the  time  evolution  operator  corresponding  to  a  Hamiltonian  (energy 
operator),  such  that  U  =  exp(-iHt),  with  U  =  cxp[-i7i(ax  ay  -  ay  ax~)/8].  Here, 
superscripts  denote  the  first  or  second  qubit  at  each  lattice  site,  and  tensor  products  of  the 
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12  1  2 

Pauli  spin  matrices  o]  are  implicit  in  the  notation,  e.g.,  ax  ay  =  ax  0  ay  .  The  collision 
operator  for  the  entire  lattice  is  a  L-fold  tensor  product  over  the  local  collision  operators, 
C  =  0j=oL"1  U. 

The  streaming  operator  S  may  be  written  as  a  product  of  operators,  each  of  which  is  a 
unitary  permutation  matrix 

10  0  0 
0  0  10 
0  10  0 
0  0  0  1 

For  a  global  shift  to  the  right  along  the  linear  lattice  of  only  the  first  qubit  at  each  lattice 
node,  S  is  a  product  over  all  but  the  right-most  site  L,  S  =  rij=iL  1  X2j-i,2j+i-  Similarly,  a 
global  shift  to  the  right  along  the  lattice  of  the  second  qubit  at  each  node  can  be  written  as 
S  =  rtj= i L' 1  X2j,2j+2-  The  inverse  of  these  matrices  are  ST,  the  transposed  operator,  and  they 
represent  a  global  shift  to  the  left  of  the  first,  respectively,  second,  qubit  at  each  lattice 
node. 

Burgers’  equation  serves  as  a  prototype  of  nonlinear  PDEs  exhibiting  front  steepening 
and  shock  formation.  As  such,  it  is  a  highly  simplified  version  of  much  more 
complicated  nonlinear  PDEs  that  have  elsewhere  been  applied  in  image  enhancement  by 
way  of  shock  front  processing  [sapiro,  osher].  The  steep  profiles  that  form  during 
Burgers’  equation  evolution  are  suggestive  of  edge  enhancement  and  image  restoration 
applications.  The  quadratically  nonlinear  term  in  Burger’s  equation  is  responsible  for  the 
front  steepening  capability. 

3.3.1  Hopf-Cole  Transformation 

Suppose  that  0(x,t)  satisfies  the  linear  heat  (diffusion)  equation  0t  =  v0xx.  If 

u(x,t)  =  -2v0x/0  =  -2v3x  In  0(x,t)  (3.23) 

then  u  satisfies  the  Burgers  equation  (B).  To  see  this,  we  consider  the  partial  derivatives 
of  u(x,t)  in  terms  of  those  of  0(x,t).  By  the  quotient  rule  we  have 

ux  =  -2v(00xx  -  0x2)/02  =  -2  0t/0  +  u2/2v,  (3 .24) 

uxx  =  -2v(00xt  -  0t0x)/02  +  uux/v,  (3.25) 

and 

ut=-2v(00xt-0t0x)/02.  (3.26) 

Then  we  have: 


(3.22) 


15 


ut  +  uux  =  -2v0xt/0  +  2v(0x/0)(0t/0)  +  uux=  -2v0xt/0  -u(0t/0)  +  uux. 


(3.27) 


We  compare  to  the  expression: 

vuxx  =  -2v0xt/0  +  2v(0x/0)(0t/0)  +  uux=  -2v0xt/0  -u(0t/0)  +  uux.  (3.28) 

Therefore  we  have  verified  that  ut  +  uux  =  vuxx,  where  the  parameter  v  is  physically  a 
scaled  kinematic  viscosity. 

We  may  note  from  equation  (3.23)  that  by  integrating,  In  0(x,t)  =  -( l/2v)Jx  u(x’,t)dx’  and 

0(x,t)  =  C(t)  exp[-(l/2v)jxx0u(x’,t)dx’],  (3.29) 

with  C(t)  =  0(xo,t).  The  latter  relation  follows  simply  from  0(xo,t)  =  C(t)  exp(O). 
Similarly  the  initial  values  of  u  and  0  are  related.  If  u(x,0)  =  Uo(x),  then 

0(x,O)  =  0o(x)  =  Co  exp [-( 1  /2v) jxx0 u0(x ’ )dx ’ ] .  (3.30) 

If  the  linear  diffusion  equation  is  solved  on  the  whole  x  axis,  then  one  way  to  represent 
0(x,t)  is  as  a  convolution  of  the  initial  condition  9o(x)  with  a  Gaussian  function  with 
variance  given  by  2vt.  Then  applying  the  Hopf-Cole  transformation  (3.23)  gives  the 
solution  u(x,t)  of  the  Burgers’  equation. 

Using  vector  notation  it  is  possible  to  write  higher  dimensional  and  vector  analogs  of  the 
Burgers’  equation  and  the  Hopf-Cole  transformation.  With  V  the  gradient  operator  in  n 
dimensions,  the  ID  Burgers’  equation  generalizes  to 

qt  +  q-Vq  =  vV2q.  (3.31) 

If  0(x,t)  satisfies  the  linear  diffusion  equation  9t  =  vV'0,  then 

q(x,t)  =  -2vVln  9(x,t)  (3.32) 

satisfies  the  n-dimensional  Burgers’  equation. 


3.3.2  Energy  Integral 

Suppose  that  we  multiply  Burgers’  equation  (B)  by  u, 

uut  +  u2ux  =  vuuxx,  (3.33) 

and  then  recognize  certain  derivatives: 
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Vi  <3tu 2  +  (1/3)  Sxu3  =  -vux2  +  vSx(uux). 


(3.34) 


Integrating  this  equation  with  respect  to  x  from  xj  to  X2  gives  the  relation 

(l/2)5t  |x2xi  u2(x,t)dx  +  ( 1/3)  [u3(x2,t)  - u3(xbt)] 

=  -lx2xi  u2x(x,t)dx  +  v[u(x2,t)ux(x2,t)  -  u(xi,t)ux(xi,t)].  (3.35) 

The  two  terms  on  the  left  side  of  this  equation  represent  the  total  rate  of  change  of  kinetic 
energy  in  the  system  and  the  net  flux  of  kinetic  energy  out  of  the  system  across  the 
boundaries.  The  two  terms  on  the  right  side  represent  the  total  dissipation  of  energy  by 
viscosity  in  the  system  and  the  rate  of  work  done  on  the  system  at  the  boundaries.  The 
nonlinear  term  of  Burgers’  equation  gives  a  means  of  putting  energy  into  the  system 
across  the  boundaries. 

We  omit  any  general  discussion  of  scaling  properties  of  the  Burgers’  equation,  but  do 
note  that  under  the  transformation  u  — »  -u  the  Burgers’  equation  becomes 


ut-uux  =  vuxx,  (3.36) 

i.e.,  the  sign  of  the  nonlinear  tenn  is  changed.  Correspondingly,  the  Hopf-Cole 
transfonnation  becomes 


u(x,t)  =  2v0x/0  =  2vdx  In  0(x,t).  (3.37) 

Since  each  term  of  the  Burgers’  equation  contains  a  differentiation,  a  linear 
transformation  of  the  dependent  variable  introduces  an  additional  linear  term  in  the  PDE. 
In  particular,  if  v(x,t)  =  au(x,t)  +  b  for  constants  a  and  b,  then  v  satisfies  the  modified 
Burgers’  equation 


vt  +  ( 1  /a)(v-b)vx  =  vvxx.  (3.38) 

This  simple  property  is  sometimes  useful  in  solving  an  initial  value  problem  by  bringing 
a  PDE  into  the  standard  form  of  Burgers’  equation. 


3.3.3  Burgers’  Equation  Simulation  Example 

We  tested  our  implementation  of  Burgers’  equation  simulation  with  a  cosinusoidal  initial 
profile 


p(x,0)  =  pa  +  pb  cos(27tx),  (3.39) 

where  pa  =  0.4  and  pb=  1.0.  The  exact  solution  for  this  problem  may  be  based  upon  a 
slightly  extended  version  of  the  Hopf-Cole  transfonnation  equation  (HC’): 

p(x,t)  =  pa  +  (2v/c)0x/0,  (3.40) 
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where  c  is  a  constant.  Then  the  solution  for  the  function  9(x,t)  may  be  developed  as  a 
Fourier  series  with  modified  Bessel  function  coefficients  depending  upon  c,  pb,  and  the 
‘viscosity’  v  [levermore].  Each  tenn  of  the  Fourier  series  expansion  exponentially 
decreases  in  time  according  to  v  and  the  square  of  the  index  of  the  Fourier  mode. 

Snapshots  of  the  dynamical  evolution  of  the  initial  profile  p(x,0)  under  the  Burgers’ 
equation  are  given  in  Figures  4a-c.  In  each  of  the  Figures  4a-c,  the  red  curve  is  the 
analytic  solution,  while  the  blue  curve  is  the  numerical  solution.  In  this  instance  the 
lattice  spacing  1  =  1  and  the  number  of  nodes  is  L  =  256.  The  very  smooth  initial  profile 
develops  a  steep  front  as  time  progresses.  The  numerical  output  of  the  QLGA  reproduces 
the  exact  solution  except  very  near  the  steep  front  where  the  solution  gradient  is  very 
large. 


Figure  4.  Analytic  and  numerical  solution  of  Burgers’  equation,  one  time  step. 


Figure  5.  Analytic  and  numerical  solution  of  Burgers’  equation,  64  time  steps. 
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Figure  6.  Analytic  and  numerical  solution  of  Burgers’  equation,  192  time  steps. 


3.3.4  Transformation  of  a  2D  Burgers'  Equation 

Very  recent  further  analysis  that  we  performed  for  Burgers’  equation  suggests  that  a 
QLGA  with  two  qubits  per  node  should  exist  for  two  and  higher  spatial  dimensions.  The 
n  =  2  case  is  again  particularly  pertinent  for  image  processing  applications.  This  section 
describes  the  candidate  2D  Burgers’  equation  we  have  identified  together  with  a  method 
for  transforming  it  to  a  known  2D  problem  whose  solution  is  known. 

Letting  subscripts  denote  partial  differentiation,  we  denote  for  example  ux  =  ou/Sx  and  dx 
=  d/dx.  We  consider  the  nonlinear  2+1 -dimensional  Burgers'  PDE 


Ut  =  -UUx  -  UUy  +  D(Uxx  +  Uyy)  +  DUx,y. 


(3.41) 


We  found  that  this  equation  should  be  efficiently  simulated  with  a  QLGA  using  the  same 
collision  operator  as  above  but  now  with  streaming  of  the  two  qubits  in  the  two 
orthogonal  x  and  y  directions.  We  have  also  found  a  transformation(s)  of  the  coordinates 
x  and  y  such  that  the  mixed  derivative  term  uxy  is  eliminated  in  favor  of  terms  ux-x'  and 
Uyy  in  the  new  coordinates. 

To  carry  this  out,  we  introduce  light  cone  coordinates 


4  =  x  -  y,  and  r\=  x  +  y. 


(3-42) 


Then  by  the  chain  rule  we  have 


3x=  ^x<9i;  +  lixSr,  =  d^  +  dn. 


(3.43a) 


8y=l)yd^  +  \\ydA  =  -d^  +  8A. 


(3.43b) 
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For  the  ’non-mixed'  second  derivatives,  by  operator  algebra  we  have 


dxx  =  (8$  +5r,)2  =  %  +  28^  +  8m  (3 .44) 

dyy  =  (-8$  +  d^f  =  %  -  2 8^  +  dm.  (3.45) 

In  particular,  for  the  mixed  derivative  term  we  have 

8x  8y  =  (8$  +  5n)(-^  +  dn)  =  8m  -  %.  (3.46) 

By  making  use  of  the  equations  just  above  in  the  2D  Burger's  equation  we  obtain 

ut  =  -2uun  +  D(u^  +  3^^).  (3.47) 


This  appears  to  be  an  anisotropic  form  of  a  2D  Burgers'  equation  very  close  to  that 
considered  by  Elton  [elton],  and  Shen  et  al.  [shen],  and  possibly  others  in  regard  to 
classical  lattice  gas  algorithms.  Using  different  lattice  spacings  Ax  ^  Ay  should 
accommodate  the  anisotropy  of  the  second  order  derivative  terms.  If  necessary,  a  scaling 
of  the  dependent  variable  u  may  be  used  to  eliminate  the  factor  of  2  in  the  2uun  tenn.  If 
u(S,ti)  =  v(£,,r|)/2,  then 


vt  =  -wn  +  D(w?  +  3v,in). 


(3.48) 


3.4  Theory  in  Support  of  the  Analysis  of  Algorithms 

Often  in  the  run-time  analysis  of  algorithms  alternating  binomial  sums  arise  in  the 
determination  of  the  order  of  operations.  In  the  course  of  our  investigations  we  were  able 
to  make  theoretical  contributions  to  this  subject  of  computer  science.  The  theory  entails 
special  functions  and  numbers  and  combinatorics  together  with  extensive  analytic 
development.  The  end  results  are  exact,  not  just  asymptotic  for  a  large  number  of 
computing  operations,  as  was  often  the  case  in  the  previous  literature.  They  support  the 
design  and  analysis  of  algorithms,  whether  they  be  classical,  purely  quantum,  or  hybrid. 

Given  an  alternating  binomial  sum,  the  sign  alternation  in  the  summand  causes 
substantial  cancellation,  masking  the  leading  behavior.  If  S(n)  =  X,  -o'1  (-1)1  ('))  f(j),  where 
(")  denotes  the  binomial  coefficient,  virtually  any  order  in  n  is  realizable,  depending  upon 
the  function  f.  If  f  =  1,  we  obtain  0,  while  if  f(j)  =  (- 1  y ,  we  obtain  2n.  Therefore,  the  full 
range  of  growth  from  sub-polynomial  in  n  to  exponential  in  n  may  result. 

Typically,  otherwise  asymptotic  results  have  been  sought  from  the  beginning  of  the 
analysis.  We  showed  how  to  exactly  calculate  a  class  of  binomial  sums  of  interest.  In 
our  approach  we  used  a  variety  of  integral  representations  and  special  numbers  and 
special  functions.  Concepts  and  functions  from  combinatorial  theory  are  additionally 
used.  Complex-valued  functions  f  are  included  in  the  theory.  For  a  recent  graduate  level 
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textbook  on  the  average  case  analysis  of  algorithms,  including  a  treatment  of  alternating 
sums,  one  may  see  the  book  of  Szpankowski  [szpankowski]. 

Our  analytic  results  are  sufficiently  general  to  subsume  many  other  previously  known 
special  cases.  From  these  more  general  results,  asymptotic  estimates  may  be  obtained  if 
desired.  Special  cases  include  a  sum  useful  for  computing  decoherence  factors  for  a 
coupled  multi-qubit  environment  and  multi-qubit  subsystem.  Our  preprint  entitled  “A  set 
of  identities  for  a  class  of  alternating  binomial  sums  arising  in  computing  applications” 
has  been  accepted  for  external  publication  [utilmath]. 


3.5  Schrodinger  Equation  Simulation 

The  time-dependent  Schrodinger  equation  (TDSE)  can  be  described  by  a  number  of 
QLGAs,  three  of  which  we  have  implemented  and  compared.  The  first  is  the  simplest, 
proposed  by  Boghosian  and  Taylor  in  their  paper  on  a  quantum  lattice-gas  model  in  d 
dimensions  [boghosian98].  The  other  two  are  improved-accuracy  algorithms  but  with  a 
greater  number  of  operations,  proposed  by  Yepez  and  Boghosian  [yepez02].  These 
QLGAs  do  not  measure  the  qubit  states  after  collision  and  streaming.  Therefore,  they 
assume  that  the  entangled  qubit  states  may  be  reliably  retained  throughout  the  length  of 
the  simulation.  In  principle  these  QLGAs  may  then  obtain  exponential  speed  up  over 
classical  simulation  of  the  Schrodinger  equation.  Because  of  the  necessity  to  maintain 
and  stream  the  entangled  qubit  states  over  all  time  steps,  the  engineering  requirements  for 
physical  implementation  are  much  more  stringent  than  for  QLGAs  with  occasional  or 
regularly  repeated  measurement.  The  Yepez  and  Boghosian  approach  uses  a  combination 
of  streaming  qubit  states  to  both  left  and  right  neighbors  in  the  case  of  a  one-dimensional 
lattice.  In  this  way  an  algorithm  more  symmetric  in  the  streaming  step  is  obtained  and 
higher  accuracy  of  the  solution  usually  results.  Versus  four  total  operations  of  streaming 
and  collision  for  the  Boghosian  and  Taylor  method,  the  Yepez  and  Boghosian  approach 
has  eight  or  sixteen  such  operations.  The  papers  cited  just  above  provide  more  details  of 
the  QLGAs. 

The  potential  energy  term  V(x)  of  the  Schrodinger  equation  is  introduced  by  applying  a 
local  phase  change  to  the  system  wave  function  as  [ibb] 

\|/(x,t)  — >  e'lV(x)T\j/(x,t).  (3.49) 

Then  we  obtain  equations  of  the  fonn 

i h  <9v|//<9t  =  (ti  2/2m)S2i|//Sx2  +  V(\j/)\j/,  (3.50) 

that  model  a  quantum  particle  under  the  influence  of  an  external  force  L(x)  =  -dV(x)/dx. 
A  key  point  is  that  we  may  take  the  potential  to  depend  upon  \j /  itself.  In  this  way,  we  are 
also  able  to  simulate  the  nonlinear  time-dependent  Schrodinger  equation.  If,  for  instance, 
V  has  a  power  law  nonlinearity  of  degree  n,  then  the  TDSE  has  nonlinearity  of  degree 
n+1.  In  particular,  taking  V(\|/)  =  |\|/|  ,  we  obtain  the  cubic  nonlinear  Schrodinger 


21 


equation.  However,  the  nonlinearities  of  the  potential  need  not  be  restricted  to  power-law 

form,  and  this  is  useful  in  modeling  other  physical  systems  [Jordan].  These  other 

2  2  2 

nonlinearities  include  V(\|/)  =  |i|/|7(l  +  |i|/|  )  and  V(\|/)  =  1-  exp(-|\|/|~). 

We  performed  TDSE  simulations  for  a  variety  of  problems  with  and  without  potential 
energy  term.  A  first  case  is  given  in  Figure  7  for  a  free  Gaussian  wave  packet.  In  this 
case  the  exact  solution  is  known  in  terms  of  separation  of  variables  and  Fourier  series  for 
the  spatial  dependence.  Another  well  known  test  case  was  provided  by  the  harmonic 
oscillator  problem.  In  Figure  8  the  corresponding  parabolic  potential  is  shown  in  green. 
The  spring  constant  k  of  the  potential  was  set  to  10°  and  the  particle  mass  to  1.  The 
classical  period  for  this  wave  packet  is  T  =  27 t/to  =  2nNk  «  1987  St,  where  St  is  eight 
time  steps  for  the  Yepez  8-step  algorithm.  The  simulation  was  run  for  80,000  time  steps, 
and  the  expected  five  periods  of  oscillation  in  the  potential  were  obtained. 


Figure  7.  Example  simulation  of  the  evolution  of  a  Gaussian  wave  packet, 

Figure  7  is  a  simulation  of  the  evolution  of  a  Gaussian  wave  packet  with  standard 
deviation  ao  =  25  on  a  ID  lattice  of  F  =  256  nodes.  The  simulated  solution  (red)  stays 
right  with  the  known  solution  (blue)  during  the  entire  simulation  over  167,000  time  steps. 


Figure  8.  Simulation  of  a  Gaussian  wave  packet  in  a  harmonic  potential. 

We  additionally  performed  simulations  of  a  Gaussian  wave  packet  encountering  positive 
and  negative  energy  barriers,  with  differing  relations  between  the  barrier  height  (or 
depth)  and  the  packet’s  energy.  For  instance,  when  the  wave  packet  energy  matches  the 
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energy  level  of  a  repulsive  barrier,  a  sort  of  resonance  occurs.  This  and  other  expected 
behavior  emerged  from  the  simulations. 

We  also  performed  simulations  for  much  more  challenging  physical  systems  including 
power-law  wave  tail  packets  (PLTPs)  and  the  Weierstrass-Mandelbrot  wave  function. 
PLTPs  have  received  much  attention  in  the  literature  in  the  last  few  years  (e.g.,  [lillo]). 

A  form  of  them  suitable  for  modeling  purposes  is 

\|/(x)  =  N  (x2  +  y2)’a/2,  (3.51) 

where  a,  y,  and  N  are  constants,  with  N  being  the  normalization  constant.  It  is  possible 
to  explicitly  write  N  as  a  function  of  a  and  y  by  evaluating  a  Beta  function  integral: 

r" 

n  = - - -  ns?'* 

^yr(a-U2W(a)’ 

where  T  is  the  Gamma  function.  We  note  that  for  large  x,  |\j/(x)|2  behaves  as  |x|‘2a. 
Hence  a  in  the  range  Vi<  a  <  3/2  causes  the  integral  (required  to  find  the  uncertainty  in 
position) 

f  x2  |  y/(x)  |2  dx  (3.53) 

J-oo 

to  diverge.  This  means  the  Heisenberg  uncertainty  relation,  while  true,  yields  no  useful 
information  about  the  uncertainty  in  momentum  or  otherwise.  In  this  case,  it  turns  out  to 
be  very  useful  then  to  investigate  entropy-like  quantities  for  the  PLTPs.  Figure  9  gives 
an  example  from  a  QLGA  simulation  of  a  PLTP  with  a  =  %  run  for  150,000  time  steps. 
This  figure  with  log-log  scale  gives  the  wave  packet  amplitude  as  a  function  of  time.  The 
packet  amplitude  decreases  with  time,  but  not  exactly  as  expected.  A  nonnalization  check 
for  the  packet  at  the  end  of  the  simulation  returned  N  =  0.966  <  1,  indicating  that  some  of 
the  packet  has  been  lost  off  the  sides  of  the  lattice  due  to  the  imposed  periodic  boundary 
condition.  This  model  problem  provided  a  stringent  test  of  the  QLGA  capabilities.  It 
presented  a  trade  off  in  the  total  number  of  lattice  nodes  to  use,  the  internal  resolution  of 
the  wave  packet,  and  the  total  number  of  time  steps  that  the  quantum  evolution  could  be 
advanced. 
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Figure  9.  Log-log  plot  of  amplitude  vs.  time  for  a  power  law  tail  packet  with  a  =  0.75. 

Finally  in  closing  this  section  and  in  regard  to  modeling  the  multidimensional 
Schrodinger  equation,  we  point  out  a  very  newly  emerging  classical  method  that  appears 
to  deserve  further  study.  This  is  based  upon  separable  representations  of  functions  in 
high  dimension  and  appears  to  overcome  the  Curse  of  Dimensionality  in  many  cases. 
The  separated  representation  pennits  many  operations  to  be  performed  with  scaling  that 
is  formally  linear  in  the  dimension.  So  we  may  identify  this  approach  as  at  least  a  good 
candidate  for  “beating  Feynman”.  This  approach  [beylkin]  now  sits  along  with  the  very 
recent  result  [dorit]  for  the  efficient  classical  simulation  of  the  Fourier  transfonn. 
Separable  representations  for  numerical  analysis  in  high  dimension  have  been  shown  to 
achieve  adequate  condition  number  and  have  been  applied  to  linear  systems  by  way  of 
iterative  methods.  In  addition,  there  seem  to  be  methods  for  dealing  with  anti-symmetric 
functions  as  arise  in  the  multi-particle  Schrodinger  equation.  A  ready  example  of  the 
latter  is  applications  in  quantum  chemistry  for  obtaining  the  ground  state  energy  of 
various  atoms  and  ions. 

To  find  a  minimal  rank  representation  for  a  separated  representation  is  still  an  open 
problem.  It  appears  to  be  a  compelling  one  for  investigation. 
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4.0  Summary 


When  it  comes  to  a  hybrid  classical-quantum  computing  environment,  quantum  lattice 
gas  algorithms  provide  a  flexible  yet  efficient  approach.  These  algorithms  combine  a 
purely  quantum  step,  the  collision  operation,  with  classical  communication  when  the 
streaming  of  qubit  states  is  performed.  Since  the  quantum  lattice  gas  algorithms  are  able 
to  model  nonlinear  systems,  they  have  a  wide  range  of  application. 

Type  II  or  hybrid  computing  offers  a  development  pathway  to  full  scale  quantum 
processors.  Hybrid  architectures  would  allow  intermediate  technology  to  be  developed 
and  demonstrated.  For  instance,  longer  coherence  times  and  error  correction  could  be 
demonstrated  within  the  array  nodes.  Besides  theoretical  and  algorithmic  work, 
experimental  effort  is  underway  on  type  II  computing  at  several  research  laboratories 
within  the  US. 

We  have  presented  algorithms  for  and  demonstrations  of  multi-dimensional  diffusion 
processing  on  a  type-II  quantum  computer  with  two  qubits  per  node.  We  presented  the 
effective  finite  difference  approximations  satisfied  by  the  density  formed  from  the  sum  of 
qubit  occupancy  probabilities  [coffeycolburn].  We  have  employed  rectangular  lattices  or 
their  higher  dimensional  analogs,  noting  that  the  lattice  spacing  in  each  Cartesian 
direction  need  not  be  the  same. 

Various  types  of  anisotropic  diffusion  have  proven  very  useful  for  image  processing, 
partly  motivating  this  investigation.  The  basic  idea  of  anisotropic  diffusion  is  to  diffuse 
intensities  along  the  edges  of  objects  that  appear  within  an  image,  while  not  diffusing  (or 
even  enhancing  the  contrast)  along  directions  that  are  perpendicular  to  edges. 

The  analysis  and  processing  of  image  data  is  an  important  and  ubiquitous  industry.  In 
addition,  diffusion,  reaction-diffusion,  and  advection-diffusion  processes  are  common 
throughout  science  and  engineering  fields. 

Of  direct  use  for  image  enhancement  purposes  is  the  ability  to  selectively  smooth  regions 
of  an  image.  We  presented  and  illustrated  an  algorithm  for  constrained  linear  diffusion  in 
two  spatial  dimensions  that  demonstrates  this  capability  on  a  type-II  computing 
architecture.  Another  direction  of  research  is  nonlinear  processing.  This  is  also  attractive 
for  QLGAs  because  they  are  able  to  simulate  nonlinear  phenomena,  whereas  purely 
quantum  methods  are  not.  Our  simulation  work  with  the  Burgers’  equation  serves  as  a 
prototype  in  this  area. 

For  two  or  higher  dimensions,  there  are  a  much  greater  number  of  diffusion  processes 
than  in  ID,  where  the  boundary  conditions  are  simply  given  at  points.  For  2D  in 
particular,  one  could  consider  other  type-II  lattices.  For  example,  it  may  be  possible  to 
employ  a  lattice  built  from  hexagons  and  with  nodes  containing  two,  three,  or  six  qubits. 
The  collision  operator  would  need  to  be  extended  to  such  situations,  and  need  not  give  a 
symmetric  mixing  of  occupancy  probabilities. 
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A  concatenated  time  sequence  of  images  provides  a  3D  space-time  in  which  to  consider 
image  processing.  By  discriminating  changing  features  in  such  data,  one  could  detect 
moving  objects.  This  has  potential  application  to  automatic  or  semi-automatic  target 
recognition. 

The  early  NMR  experiment  which  qualitatively  demonstrates  the  feasibility  of 
implementing  ID  diffusion  [pravia]  is  encouraging.  It  provides  evidence  that  the  path  to 
large  hybrid  quantum-classical  computers  may  be  realizable  well  before  large-scale  type- 
I  quantum  computers.  In  addition,  liquid  state  NMR  realizations  of  an  ensemble  of 
sufficient  size  would  exceed  the  computational  capacity  of  classical  digital  computers. 
The  various  diffusion  attributes  of  dimension,  isotropy  or  anisotropy,  forward  or 
backward,  and  linear  or  nonlinear  processes  provide  opportunities  for  designing  aspects 
of  image  processing  on  a  type-II  quantum  computing  architecture. 
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